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ABSTRACT 

In order to interpret the results of complex realistic star cluster simulations, which rely on 
many simplifying approximations and assumptions, it is essential to study the behavior of even 
I more idealized models, which can highlight the essential physical effects and are amenable 

to more exact methods. With this aim, we present the results of iV-body calculations of the 
evolution of equal-mass models, starting with primordial binary fractions of - 100 %, with 
values of N ranging from 256 to 16384. This allows us to extrapolate the main features of the 
^— ^ ■ evolution to systems comparable in particle number with globular clusters. 

\^ I In this range, we find that the steady-state 'deuterium main sequence' is characterized 

. by a ratio of the core radius to half-mass radius that follows qualitatively the analytical esti- 

mate by Vesperini & Chernoff ( 1994), although the dependence is steeper than expected. 
Interestingly, for an initial binary fraction / greater than 10%, the binary heating in the core 
during the post collapse phase almost saturates (becoming nearly independent of /), and so 
O ' little variation in the structural properties is observed. Thus, although we observe a signifi- 

^ . cantly lower binary abundance in the core with respect to the Fokker-Planck simulations by 

^2 ' lOao et al. ( 1991 ), this is of little dynamical consequence. 



OO 
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At variance with the study of Gao et al. (1991), we see no sign of gravothermal oscil- 
lations before 150 halfmass relaxation times. At later times, however, oscillations become 
prominent. We demonstrate the gravothermal nature of these oscillations. 



5-H ' Key words: globular clusters - methods: N-body simulations - stellar dynamics. 



1 INTRODUCTION 

The evidence for primordial binary stars in old globular star clus- 
ters continues to accumulate st eadily. Even thoug h no comprehen- 
sive re view has appeared s ince |HuLeLalJ (l99^. recent additions 
includ e iPulone et alJ ^2003^ ■ lBellazzini et alJ <2002h . lAlbrow et alJ 
j200lh . 

The A''-body modeling of such clusters has also devel- 
oped rapidly, with studies devoted to exotic species associ- 
ated with binaries, blue stragglers, intermediate mass black 
holes, hierarchical systems, and so on (see, for example, 
ftlurleyetal. 2001; Portegies Zwart & McMillan 2002; Aarseth 
BoOltlPortegies Zwart et a l. 2004; H urley et alj2005lh . 

There have been relatively few A'^-body studies devoted to the 
effects of binar i es on t he structural evolution o f model clusters. 
[McMillan et al.' ('199(f), 'McMillan et al." ('l99l'), 'McMillan & Hut' 
a9941 . .Heggie & Aarseth I.19921 . .Aarseth & Heggie. L1993.) . 
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iKroupd 11995^ and'de La Fuente Marcos! |l99?) described a wide 
range of results with values of N which are generally rather 
modest by modern standards (though not inappropriate for appli- 
cation to many open clusters). More recent studies include that 
of Wilkinson et al. (2003), whose focus was on the evolution of 
the core radius with age. Generally speaking, however, in recent 
years the effect of primordial binaries on the structural evolu- 
tion of rich systems has been studied using approximate methods 
based on the Fokker-Planck equation or related approaches. The 
study by Gao et al. ( 1991) [hereafter Gao et al.], using the Fokker- 
Planck model with finite differences, establis hed an interesting 
set of parameters which has been taken up bv lGiersz & SpurzemI 
(2000, 200^, using their hybrid gas/Monte Carlo code, and by 
iFregea u et al. (2003, 2005), whose code was pure Monte Carlo 
(with computed binary interactions in the later paper). 

For all the usual reasons, it is advisable to check the results 
of approximate methods by means of A'^-body simulations, and it 
is the main purpose of the present paper to do so, based on the 
standard model defined by Gao et al. We present the results of an 
extensive set of direct A'^-body simulations, performed with up to 
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16384 particles and different ratios for the number of initial bina- 
ries to single stars, ranging from 0% to 100%. Following an outline 
of the physical picture in Sec.|2| the specification of our simulations 
is given in Sec.|3| 

The subsequent two sections contain the main presentation of 
our results, with emphasis on (i) a comparison with the results of 
Gao et al. (Sec.|4j, and (ii) an analysis in terms of the dependence 
on A'^ and on the primordial binary fraction (Sec.|5}- In fact these 
topics cannot be cleanly separated. For example, the main system- 
atic differences observed with respect to Gao et al. are in the frac- 
tion of binaries in the core, and possibly in the size of the core, 
in the post-collapse phase. After extrapolating our results to the 
number of particles used by Gao et al., we clearly observe a lower 
fraction of binaries in the core (by around a factor 3). Extrapolation 
of the core radius is more problematic, however (see Sec. 15.1.^ . 
In this section we also compare our results with those of other ap- 
proximate numerical methods, and with the analytical estimate of 
the core radius due to lVesperini & C hernofi ( 1994). The last sec- 
tion of the paper provides an extended summary. 



2 THE PHYSICAL PICTURE 

For an isolated cluster, the most important characteristics describ- 
ing the mass distribution are the core radius rc and the half-mass 
radius ru. Fig.Qshows how both evolve for the case of a Plummer 
model distribution of single stars, all of which have equal mass, 
with no primordial binaries present. The half-mass radius remains 
roughly constant during the core collapse phase, while the core 
shrinks at an increasing rate until the first core bounce takes place. 
Around this time, one or more binaries are formed dynamically, 
through simultaneous three-body encounters. The energy produc- 
tion in those binaries halts core collapse, and also initiates an ex- 
pansion of the half-mass radius. In the post-collapse phase, the core 
expands and contracts dramatically, while new binaries are formed 
and older binaries are first hardened and then ejected. 

Fig. |2| shows the striking difference that primordial binaries 
make to this simple picture. The presence of primordial binaries 
causes the core to shrink roughly twice as fast as in the single-star 
case. This effect is largely due to mass segregation: the binaries, 
between twice as massive as single stars, tend to precipitate toward 
the center. In contrast, the core contraction is much less deep in the 
case of primordial binaries: the core contraction levels off after the 
core shrinks by less than a factor ten, whereas in the single star case 
core collapse continues till the core is one hundredth of its original 
size. Another difference is that in the primordial binary case, the 
half-mass radius starts increasing very early on, a factor four 
earlier than in the single star case. Since there is very little mass 
loss in these early stages of cluster evolution, an increase in half- 
mass radius indicates an injection of energy into the bulk motion of 
the stars and the center of masses of the binaries. This points to an 
efficient energy mechanism at work, which then also explains the 
early halting of the core contraction. 

This is indeed the physical explanation. Binaries can give off 
energy in encounters between binaries and single stars^ (and also in 
binary-binary encounters), the probability of which is given by the 
product of the densities pbps of binaries and single stars, whereas 

^ Though such encounters can result in the ejection of single stars and bi- 
naries, the cluster still gains energy, by the process known as "indirect heat- 
ing". 
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Figure 1. Dependence of half-mass (upper curve) radius (N-body units) and 
of the ratio rc/r^ on time (units of the initial half-mass relaxation time). 
The simulation has been performed with 16384 particles and no primordial 
binaries. 

the formation rate of new binaries is proportional to . The rel- 
evant scales for these densities are set by the ninety degree turn- 
around distance d, which is of order 1/Af in standard units, far 
smaller than typical inter-particle distances, even in the core. If we 
define the interaction densities as the number of particles per unit of 
volume dp , and denote them by a hat, we have /5s ^ 1 and p(, ^ 1. 
In a model with primordial binaries, ps and p(, are comparable, at 
least in order of magnitude, and it follows that phPs S> p^. In other 
words, at a given density the rate of energy generation by primor- 
dial binaries far exceeds the rate of energy generatio n by binaries 
that first have to be formed at that density I Hut 198?). 

Another way of describing the difference between these two 
forms of evolution is that Fig. Q shows a collapse to the 'hydrogen 
main sequence', whereas Fig.|2|shows a collapse to the 'deuterium 
main sequence', to borrow metaphors from stellar evolution. Let 
us take a gas cloud of a solar mass, consisting of only hydrogen 
and helium. When this protostar contracts sufficiently, nuclear re- 
actions in the core will at some point balance the radiation losses 
at the photosphere. At this point, by definition the star has landed 
on the (hydrogen) main sequence. If the star would be made of 
pure helium, it would shrink further, until it would reach the he- 
lium main sequence, where the central temperature and pressure 
are high enough to generate the same energy at the rate at which it 
is lost at the surface. In practice, however, protostars never contain 
pure hydrogen, but typical have a very small admixture of deu- 
terium. This deuterium can bum to form helium at lower densities 
and temperatures than hydrogen can, for a similar reason that pri- 
mordial binaries can 'burn' gravitationally at lower densities than 
single stars. As a result, a contracting protostar first gets hung up 
at the deuterium main sequence, where it has a radius significantly 
larger than it will have soon afterwards, when its central deuterium 
reservoir is burned up, and the star settles on the (hydrogen) main 
sequence. 

Here we witness something similar: our contracting cluster 
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Figure 2. Dependence of half-mass (upper curve) radius (N-body units) and of tlie ratio rc/r;, on time (units of the initial half-mass relaxation time). The 
simulation has been performed with 16384 particles and 10% binaries (as defined in Eg. IT]. 



core settles onto an almost stable extended core configuration, 
while burning the primordial binaries. Gradually, when this primor- 
dial fuel is being exhausted, the core shrinks since the burning pro- 
cess becomes less efficient when the binary fraction becomes less. 
This slow shrinking of the core continues until (almost) all primor- 
dial binaries are burned up, after which the new binaries have to be 
formed by the process of simultaneous three-body encounters be- 
tween single stars. Note that the Tc value in Fig.|2|toward the end is 
dropping toward the values seen in Fig.Qfrom the moment of the 
first core collapse. This behavior is more evident in Fig|5| 



3 SIMULATIONS: SETUP AND ANALYSIS 
3.1 Initial conditions and definitions 

The models shown in Figs. Q and |2| like all those studied in this 
paper, are isolated, with stars of equal mass. The initial distribution 
is a Plummer model. In the standard model of Gao et al., initially 
10% of the objects are binaries, as in Fig.|2| with energies in the 
range 3 fcTc-400 fcTc, where Tc is the central temperature of the 
system. For a Plummer model Tc ~ 1-7 T, where (3/2)fcT is the 
mean kinetic energy per particle of the system (the binaries being 
replaced by their barycenter). Gao et al. used Fokker-Planck models 
with parameters corresponding to A'^ = 3 10^ particles to study the 
evolution. They also carried out runs with the initial fraction of 
binaries equal to 0.05 and 0.2. 

In this paper we report on the result of A'^-body calc ulations, 
which we performed using Aarseth's NBODY6 code iAarsethI 
l2003h . which has been slightly modified to provide additional di- 
agnostics on the spatial distributions of single and binary stars (see, 
e.g.. Fig, not . Starting from a minimum of 256 particles, we have 
made runs with no primordial binaries, 2% primordial binaries, the 
three Gao et al. choices of 5%, 10%, and 20% primordial binaries, 
as well as runs with 50% and 100% primordial binaries (Table0. 
As in the case of Gao et al., we have defined here the primordial 
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Figure 3. Evolution of the core radius versus time for different amounts 
of smoothing. From top to bottom the length of the smoothing window is 
10trh(0),ltj.h(0), Itii, lO^^td- For all but the last one, the curves are 
shifted upward by a constant factor for a better graphical representation. 
The simulation has been performed with 1024 particles and 10% binaries. 



binary fraction as the fraction of objects that are binaries: 

f ^ Nb/{Ns+Nt,) 



(1) 



where and Ns are the initial number of binaries and single stars, 
respectively. This implies that the fraction of the total mass in bi- 
nary stars, in the case of equal masses, is larger in the following 
way: 

= 2Ni,/{N, + 2iVf,) 
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Table 1. N-body models with Gao et al.-like initial conditions 



N f: 


0% 


2% 


5% 


10% 


20% 


50% 


100% 


256 


43 


2 


2 


46 


39 


16 


2 


512 


49 


2 


3 


46 


45 


24 


1 


1024 


51 


3 


4 


45 


44 


46 


1 


2048 


47 


3 


2 


36 


31 


21 


2 


4096 


9 


3 


2 


13 


8 


8 


2 


8192 


1 


1 


1 


3 


1 






16384 


1 






1 









Note: This table gives the number of realisations for each value of TV and 
/■ 



For example, for a run with 10% primordial binaries, / — 0.1 
whereas = 2/11 ~ 0.18. 

Note that A'^ denotes the number of original objects, i.e. A*' = 
Ns + Nb. When we discuss a run with = 1024 and 50% pri- 
mordial binaries, we are dealing with 1536 stars. 

All our results are presented using standard units 
jHeggie & Mathietll98d) in which 



G = M 



-4Et 



where G is the gravitational constant, AI the total mass, and 
Et the total energy of the system of bound objects. In other 
words, Et does not include the internal binding energy of the 
binaries, only the kinetic energy of their center-of-mass motion 
and the potential energy contribution where each binary is con- 
sidered to be a point mass. The corresponding unit of time is 
td = (GM)^^^ /{-4:Et)^^^ = 1. For the relaxation time, we use 
the following expression 



trh ~ 



0.138iVr 



3/2 



In(O.lliV) 

A quantity of great interest is the core radius, which we here 
define as in NBODY6, i.e. essentially 



X]i=l,JV Pi 

where is the distance of the ith star from the density centre, and 
the density pi around each pa rticle is computed from the distance 
to the fifth nearest neighbour iCasertano & Hullll985l) . In practice 
outlying stars are omitted from the sum. 

The runs were performed using idle time on many worksta- 
tions through the Condor system^: this allowed us to obtain in a 
few months what would have taken a few years on a single dedi- 
cated workstation. 



3.2 Fluctuations 

One significant difference between the A'^-body models we are us- 
ing and the models of Gao et al is the presence, in A'-body models, 
of fluctuations. These have both a cosmetic and a scientific impor- 
tance, and we consider these issues together here. 

In presenting the results from A/^-body simulations, there is 
always the question of how much smoothing to apply. In Fig. 3, we 
present the results for one run, for the change in time of the core 
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Figure 4. Change in internal energy of bound binaries, per unit time. The 
results have been smoothed with two triangular windows of full width 
2 t^h(O) (solid) and 10 t^^{0) (dotted). 



radius, for a variety of different amounts of smoothing. In some 
figures in this paper we apply smoothing in order to highlight the 
particular points made in each figure. Fig.|3|can be used as a gauge 
by which to judge the amount of 'noise' present in various figures. 
Generally we have used a smoothing interval of 2.5ir-h(0). 

These fluctuations also have a scientific interest. In the con- 
text of three-body binaries (i.e. those formed in three-body en- 
counters, without a population of pri mordial binaries ) their role 
was considered in simi jlified models bv llnagaki & Hull jl988h and 
Ijakahashi & Inagaki' ('1991'). The role of fluctuations in an A^- 
body rnodel, though in the absence of binaries, was considered by 
Isettwieser & SugimotoHl983) . They showed that values of the ve- 
locity dispersion, sufficiently well defined to determine its spatial 
gradient, could be determined in a 1000-body system, if averaged 
over a shell containing at least 90 particles and over at least one 
crossing time. With our present models we have an opportunity to 
consider similar issues for the generation of energy by primordial 
binaries. 

We restricted attention to our largest simulation (TV — 
16K, f — 0.1) and considered the change in internal energy of 
all the bound binaries. Within one A'-body time unit this is usually 
zero, but may be negative (if a binary hardens, for example) or pos- 
itive (if a binary escapes). Perhaps surprisingly, it is necessary to 
smooth over a long time interval (at least 2trh{0)) in order for the 
result to have a consistent sign (Fig|4}. 

We have also considered how these fluctuations affect the rate 
of interactions, in the following way. If the proportion of binaries 
in the core is very high, the rate of binary-binary interactions is 
Nbb oc PcTc, where pc is the central density. Since pcV^ oc 
it follows that Ntb oc vi/rc- Therefore a simplified model, such 
as a Fokker-Planck model, in which binaries interact at the cor- 
rect average rate, will yield a core with parameters vfpjVfp such 

that = (— ■). Fluctuations in an A^-body model, however, will 

rpp rc 

lead to values of the mean core radius and central velocity disper- 
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sion such that 



and so one need not expect these mean values to agree with those 
of a Fokker-Planck model, even if the rate of energy generation by 
binaries in the Fokker-Planck model is correct. However, we have 
checked that the difference in eq.j2j, with smoothing as in Fig|4| is 
very small. In a core supported by primordial binaries, the core con- 
tains many particles, and so fluctuations are unimportant by com- 
parison, say, with a core supported only by three-body binaries. 



4 COMPARISON WITH GAO ET AL. 

The main purpose of this section is to present results of A'^-body 
simulations which may be compared rather directly with those in 
Gao et al. Therefore the sequence of presentation closely follows 
the figures in that paper. The main difference between our runs and 
those of Gao et al. is in the particle number. While Gao et al. used 
parameters appropriate to A'^ = 3 x 10^, even our largest run is 
considerably smaller. Our discussion of the Af-dependence of our 
results is mainly given later, in Sec.|5| T he runs discussed by Gao e t 
al. have also been taken as a test case by ' Giersz &Spurzenjt2OO0l) 
and bv Ipreg eau et al. ( 2003, 2005), and we also make occasional 
brief comparisons with their results. 



4.1 Core and Half-Mass Radius 

The structural evolution is best captured by the half-mass and core 
radii. Fig.|2|shows the results for a run with A'^ = 8192 and 10% 
primordial binaries, and can be compared with Fig. 1 in Gao et al., 
except that they scale radii by the initial scale radius of the Plum- 
mer model (i.e. about 0.59 in our A'^-body units). The expansion 
of Th is similar (about 0.7 dex in both calculations), and the main 
differences are in the evolution of the core radius. These are 

(i) in the A'^-body model the initial decrease of rc is smaller, and 
T'c remains considerably larger than in the Fokker-Planck model; 
and 

(ii) in the A'-body models (even for A'^ = 16384), there is no 
sudden onset of deep core collapses, as observed in the Fokker- 
Planck model at about t/trh {0) — 50. 

We now discuss these two issues in turn. 

The fact that the contraction of the core is less deep than in the 
model of Gao et al. could partly be an effect of the A'^-dependence 
which we observe in our A'-body models (Fig. llSt . Our conclusion, 
however, depends on how the extrapolation to A' = 3 x 10^ (the 
va lue used by Gao et al.) is ca rried out. If we adopted the theory 
of IVesperini & Chernofdil994l (see Equation |3} we would expect 
Vc/rh — 0.04, which still exceeds by about 60% the value of about 
0.025 found by Gao et al. The data from our simulations, however, 
suggest a steeper decrease of the core size with A', and so the num- 
ber given by Gao et al. could well be consistent with our direct 
simulations. 

The difference in values of could be caused by other fac- 
tors, including the relative abundance of single and binary stars 
(Fig. II H . which is also A'^-dependent (Fig. I21> . A further factor 
to be borne in mind is that the core radius should depend on the ef- 
ficiency at which binaries produce energy in three- and four-body 
collisions, and Gao et al. were well aware of the uncertainties in 



Table 2. Comparative results on the core radius and gravothermal oscilla- 
tions for simulations with / = 10 % 



Source 


i5 logio r-c 


tgto 


This paper 


0.65 


130 


^iersz & Snurzem (2000) 


0.2 


160 


Giersz & Snurzem (2003) 


0.42 


20 


Gao et al. (1991) 


1 


50 


Freaeau et al. (2003) 


0.48 


70 


Freaeau et al. (2005) 


1 


^ 125 



Notes: S logj^Q r,, is the change in the logarithm of rc from t = until the 
end of the initial contraction onto the "deuterium main sequence"; tgto is 
the time (in units of t,.h(0)) of onset of gravotheiTnal oscillations. Most of 
these values have been estimated quite approximately from the figures in 
the quoted papers. 

the analytical cross sections which they adopted. Finally, the defi- 
nition of core radius in the two types of simulation may be different, 
though the initial values are quite similar. 

Comparison with other results in the literature is beset by simi- 
lar uncertainties. We report simply that Fregeau et al. 1 2003) found 
an initial contraction of the core fairly co mparable to our resul t 
(despite the different values of N), but th enlFreeeau et all i2005h 
obtained a result comparable with that of iGaoe^LTl^^l) when 
binary interactions were computed directly. In the hybrid model of 
Foiers z & Smirzem 1 2000), on the other hand, the initial decrease of 
the core radius was only about 0.2 dex (Table|2j. 

Now we consider the onset of gravothermal oscillations in 
our simulations. While gravothermal oscillations would not be 
expected to occur in models as small as that shown in Fig. |2| 
(Xjoodman 1987), for A^ — 16384 deep gravothermal os cillations 
are o bserved in single-component A'^-body systems (cf. iMakind 
Il99d) . Indeed we observe gravothermal oscillations in our run with 
A'^ = 16384. However their onset occurs later than in the simula- 
tions of Gao et al. (see Figs.|5||6j, though again we do not know 
how the time of onset should vary with A'^; this difference may or 
may not be significant. 

Further evidence, though it is contradictory, is provided by the 
resul ts of other authors ( Table O. By comparison of their two pa- 
pers, [oiMS^^^purzem' (2003) concluded that the time of onset 
of gravothermal oscillations is greatly affected by the treatm ent of 
binary interactions. The treatment of iFregeau et alj i2003h more 
closely resembles that used by Giersz & Spurzem in their ear- 
lier paper, but they found (Fregeau et al. (2005)) that the onset of 
gravothermal oscillations is delayed when binary interactions are 
computed directly. 

4.2 Total Mass in Singles and Binaries 

The destruction of the binary population is shown in Fig.Q which 
is comparable with Fig. 2a in Gao et al. The destruction of binaries 
takes place at a quite similar rate in the two models, but escape of 
single stars is much larger in the A'-body simulations. Indeed in 
the A'^-body models this somewhat diminishes the early increase in 
the numbers of single stars which is so prominent in the Fokker- 
Planck model, and which is caused by the destruction of binaries; 
nevertheless this f eature is still visible in our A'^-body results. In 
the hybrid model of lGiersz & Spurzem ( 2000) the number of single 
stars decreases monotonically, though the residual mass in binaries 
at t = lOOtrh (0) is very similar to the value in our Fig.Q 

The cumulative mass lost by escape (separately in binary and 
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Figure 5. Gravothermal oscillation in the runs with A'^ = 16384 and 
10% primordial binaries run (continuous curve) and no primordial bina- 
ries (dashed line). The data have been smoothed with a smoothing interval 
2.5 UiM- 



Figure 7. Time-dependence of the number of single (upper dotted line) and 
binary stars (lower dotted line) expressed as a fraction of the initial values. 
The solid line is a graph of the total mass of the system. The data refers to a 
simulation with = 16384 and 10% primordial binaries. The unit of time 
is the initial half-mass relaxation time. 




Figure 6. Three successive clockwise gravothermal oscillations in the N = 
16384 10% primordial binaries run. The number of particles in the core 
is plotted against the core radius. An isothermal core is a two-pai'ameter 
family, and the presence of loops in such a diagram may be regarded as 
diagnostic of oscillations driven by gravothermal effects (Makinc 1996i). 
The three main loops, which are distinguished by different line types for 
clarity, correspond to the last three complete oscillations in Fig.lsl 




ao 40 60 80 
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Figure 8. Evolution of the escaper mass in singles (Mes, solid line) and 
binaries (Mgi,, dotted line) for a simulation with A'^ = 16384 and 10% 
primordial binaries. 

single stars) is shown in Fig.js] which can be compared with Fig. 27 
in Giersz & Sourzem (2000). The mass of single escapers at time 
lOOtrh(O) is well matched, but we find that the mass of binary 
escapers is about 50% larger than in the hybrid model. 

In our model with A'^ = 16384, the decrease of the total mass 
by time t — lOOtrh (0) is rather similar to that found by Fregeau et 
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al (2003, Fig.4; 2005, Fig.2), but is about 50% larger in our model 
than theirs by time t = 200trh(0). We find that the total potential 
and kinetic energy of the bound members varies with t in a man- 
ner very similar quan titatively to the result depicted in their Fig. 2 
jFregeauetai]j2003h ^. 



4.3 Energy Distribution of Binaries 

Fig. |9| shows the evolution of the distribution of internal binding 
energies of the binaries, and can be compared with Fig. 2b in Gao 
et al. Both diagrams use a logarithmic energy scale (abscissa), but 
with different definitions; for Gao et al. it is based on the discrete 
energy levels they adopted, whereas in our case it is determined by 
the binning in the output of NBODY6. The range of energies shown 
in the figure of Gao et al. (which is slightly smaller than the range 
of energies in the initial conditions) corresponds in our diagram to 
the range 2.1 ^ Iog2(-Ei,/-Eso) ^ 8.5, where Ebo is defined in 
the caption. 

Harder binaries were not plotted by Gao et al., as they were 
treated as being dynamically inert. By the end of their run this 
group of binaries accounted for about one third of the remaining 
binaries. In our A'^-body simulations, very few such extremely hard 
binaries are found. Simple theoretical ideas show that a binary is li- 
able to be ejected (in a binary-singJe interaction) if Ei, iZj 15m\(f>c\ 
jGoodmantl984^ . where cpc is the central potential. For a Plummer 
model this gives Ei,/Ebo ~ 100, which corresponds to a value of 
about 7 in the units of the horizontal axis of Fig.|9| in good agree- 
ment with our results. The central potential deepens by a factor of 
only about 2 (Fig. ll5> . 

We consider that the escape of very hard binaries in the A'^- 
body simulations accounts for the main differences between our 
result (Fig.|3 and that of Gao et al. While we find that the distribu- 
tion of Et (at fixed time) declines at the highest binding energies, in 
the results of Gao et al. the distribution of binding energy is mono- 
tonically increasing. Both kinds of simulation agree in the deple- 
tion of less-hard pairs, though in Gao et al. the depletion extends to 
somewhat harder pairs. To quantify these effects briefly, we remark 
that, for A'^ — 16384, the lower 10 percentile increases by a factor 
of about 3.5, and the mean increases by a factor of about 2.3, which 
is smaller because of the preferential destruction of softer binaries. 

4.4 Spatial Distribution of Binaries 

The evolution of the spatial distribution of binaries is indicated in 
two ways in Figs. llOl and ll II which can be compared with Fig. 3 in 
Gao et al. 

Our results for the evolution of the half-mass radius of the sin- 
gle an d binary stars (Fig. |10t are very similar to those of Gao et 
al. and lOiersz & Spurzem i2000), and rath er at variance (at leas t 
in later stages of the evolution) with those o f lFregeauetJ]j2003h : 
their Fig. 4 indicates the half-mass radius of the binaries exceed- 
ing its initial value by about 70trh{0), and exceeding the cur- 
rent half-mass radius of the single stars after about 120t,.fe(0). In 
iFregeau et alJi200S') the evolution of the half-mass radius is, if any- 
thing, still faster. 

Within the core the agreement with Gao et al. is poorer. Our 
values for the density ratio (in binaries and sing les: Fig.fm are con- 
siderably smaller than the values found by Gao et al.; those authors 
find po,b/po,s — 10 towards the end of the steady binary burning 
phase. This quantity is A^-dependent (see Sec. 15. 2. 2t . but we do not 
consider this sufficient to explain the disagreement. (Extrapolation 
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Figure 9. Evolution of the distribution of binary binding energy (given as 
in NBODY6 in units of Ebo = {3/2)kT, where T is the mean kinetic 
temperature of the initial Plummer model). The ordinate is the number of 
binaries per bin. which are of unit length in log2- Curves have unit normal- 
ization, and are shown for 4/4^^(0) = 0, 20, 40, 60, 80. The simulation 
has been performed with 8192 particles and 10% of binaries. 
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Figure 10. Dependence of half-mass radius (N-body units) for singles (up- 
per curve), whole system (second curve) and binaiies (central curve) on 
time (units of the initial half-mass relaxation time). The lowest curve is the 
core radius of the system. The simulation has been performed with 8192 
particles and 10% of binaries. 
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Figure 11. Binary fraction in tlie core, as a function oft/t^i^ (0). Tlie upper 
curve is the ratio of tlie mass density of binaries to the mass density of single 
stars; the lower curve is the binary fraction (corresponding to /, eg.fTl. 
but restricted to the core). The simulation has been performed with 8192 
particles and 10% of binaries. 



from our data suggests po,b/po.s ~ 3 for the Gao at al. number of 
particles). 

From the dynamical point of view the early core contraction 
and increase in pt / Ps can be seen as manifestations of mass segre- 
gation. Soitzet 1 1969) showed that this could lead to equipartition 
between two species (here binaries and single stars) if the initial 
mass in the heavier species was sufficiently small. His criterion cor- 
responds, in the present situation, to / 0.03. Therefore all of the 
evolution we observe takes place in the regime corresponding to 
the mass segregation instability. 

It is worth considering whether the topics discussed in this 
subsection and the last are interconnected, i.e. whether the distribu- 
tion of binding en e rgies is different at different radii. Fig. 24 of 
iGiersz & S ourzem (2000) implies that the removal of soft pairs 
mainly affects binaries in the core. They also find a gap in the 
numbers of bina ries at intermediate radii, which was predicted by 
|Hutet^lfT99?). These trends seem to be confirmed in our simula- 
tions (see Fig. ll2l bottom right panel), but the number of surviving 
binaries is too low to draw firm conclusions. 



4.5 Lagrangian Radii 

We have already discussed the evolution of the half-mass radius. 
The evolution of other Lagrangian radii is displayed in Fig. 1131 
which can be compared with Fig. 5 in Gao et al., except that the 
choice of mass fractions differs at some points from theirs. The 
differences in the Lagrangian radii are at the level of quantitative 
detail. For example, we find that the 10% Lagrangian radius grows 
by about 0.1 dex up to f — 50trh(0), whereas Gao et al. find an 
increase by about 0.3 dex. 
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Figure 12. Evolution of the distribution of binaries in space of radius 
and binding energy. The radius is given with respect to the instantaneous 
core radius; the energy in units of kT. The panels refer to t/t,.i^{0) = 
0, 15, 45, 90. The dashed line is the position of the half mass radius. The 
simulation has been perfor med with 4096 particles and 10% binaries. Com- 
parison with similar data in iGiersz & SpurzenJ l200dl) is given in the text. 
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Figure 13. Lagrangian radii (enclosing 2, 10, 20, 50, 75, 85% of the total 
mass), for the simulation with A'^ = 16384 and 10% of primordial binaries. 

4.6 Structural Parameters 

Gao et al. point out that much of the evolution of their model can 
be described approximately as homological. Part of the evidence 
for this is their finding that some dimensionless measures of their 
model vary little over the evolution, even though individual quanti- 
ties, such as the half-mass radius, may change by large factors. We 
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Figure 14. Structural parameter for the star clusters ETr^/GM^, where 
i?x is defined in Sec.|3] To be compared with Fig 6 in Gao et al (though they 
plot separate curves for AIcpo/E-Y and t^E-y / {GKP)). The simulation 
has been performed with 8192 particles and 10% of binaries. 



Figure 15. Structural parameter for the star clusters tfuWh/GM . To be com- 
pared with Fig 6 in Gao et al. The simulation has been performed with 8192 
particles and 10% of binaries. 



find, as do Gao et al., that the structural parameter E^ytu/GM^ is 
almost constant (Fig. 1141 note the magnified vertical scale). 

The value of —0.2 is often adopted for this scaled energy in 
simple models of bound stellar systems ( Spitzgt 1987. equation (1- 
10)). 

Another important quantity is the central potential, a scaled 
form of which is presented in Fig. 1151 Though our scaling differs 
from that adopted by Gao et al., the fractional change between the 
initial value and the value in the post-collapse phase is comparable. 



5 DEPENDENCE ON N AND BINARY FRACTION 
5.1 Core evolution 

5.1.1 The time of core collapse 

Determination of the time of core collapse requires first a care- 
ful operational definition of this term, which is much less clear in 
these models than it is in the absence of primordial binaries. If A'^ is 
large enough, when / = there is a deep and sharp core collapse 
(Fig.0, and the determination of the time of core collapse is rela- 
tively un ambiguous. It does, howev er, vary from one simulation to 
another jSpurzem & Aarsethlll99d) . In the presence of primordial 
binaries the problem is illustrated in Fig. ll6l The initial collapse of 
the core is almost linear, and we may determine the point at which 
this ends {tec) by finding the intersection of two linear fits to the 
data, for t > tec and t < tec respectively. This has been done itera- 
tively. Alternatively, one may attempt to determine the end of core 
collapse by finding the minimum of rc. This requires smoothing 
of the data, and in general the result will depend on the smoothing 
interval. 

Both methods present advantages and disadvantages. The lin- 
ear fitting method is very robust for the determination of the first 
line (i.e. for t < tec), which plays the major role for setting the 



time of core collapse. In fact this first line is much steeper than the 
second (that fits the core radius in the phase of post collapse). This 
second fit is more sensitive to the noise and to the extent of the 
fitting window (and thus in particular to the end time of the simu- 
lation), so that the uncertainty associated with the determination of 
the core radius at the time of core collapse, defined as the intersec- 
tion of the two linear fits, is larger than the uncertainty in the time. 
It may be also argued that, with this definition, the core radius de- 
termination is somewhat artificial, since it does not correspond to 
a realized (or averaged) value in the simulation. In fact, at the time 
of core collapse the core radius in the simulation is bigger than the 
intersection of the two fits (see Fig. ll6> . 

The second method, i.e. the determination of the time of core 
collapse as the time at which the minimum value of rc is attained, 
may be significantly affected by the intrinsic random fluctuations 
in the measure, due to its sensitivity to the local details of the rc{t) 
curve. After some experimentations, we found that a good choice 
for minimizing the errors is given by adopting a smoothing win- 
dow with a width of 2 - 3 trh (0) for iV > 2048 and slightly bigger 
(3-4 trh(O)) for iV < 2048. We note, however, that the local 
nature of this estimate can play also a welcome role, since to de- 
termine the time of core collapse, and the associated core radius, 
the simulations may be stopped just shortly after core collapse. In 
contrast, the double linear method requires data for several tens of 
relaxation times during the postcollapse phase. 

To evaluate the relative performance of the two estimators we 
have applied both methods to a large number of simulations (a ma- 
jor subsample of the simulations with 20% primordial binaries; 
see Table 1). As can be guessed from Fig. 1161 we found a sys- 
tematic bias of approximately 3trh(0) between the two methods 
in the measure of the time of core collapse, with the linear fit giv- 
ing smaller times. The numerical experiments also confirm that the 
variance of the measure is marginally smaller (30% less) for the lin- 
ear fit method; the results are only very weakly N-dependent. The 
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Figure 16. Determination of tlie time of core collapse. Smoothed and un- 
smoothed core radius for the model with A'^ = 2048 and 10% binaries. The 
two linear graphs are explained in the text and are used to determine the 
time for core collapse. 



values for the core radius at core collapse are, on the other side, 
similar for the two methods. 

In what follows we present results from the second method 
(minimum of Vc) for determining t^c and r^c, and the possible lim- 
itations of its significance need to be borne in mind. The main 
reason for choosing this procedure is the fact that this has al- 
lowed us to consider as valid simulations all those that reached at 
least 10 trh(O) after core collapse. In fact, while formally all the 
runs should have lasted for at least lOOtrh(O), some of them were 
stopped prematurely due to numerical difficulties. Given the large 
number of experiments (see Table 1), it was impractical for us to 
manually restart all the runs that required assistance. For each en- 
try in Table 1 we have, however, at least one completed simulation, 
and on average the completeness of the sample is above 75%. 

Our results are shown in Fig. 1171 From this figure it may be 
concluded that, at a given value of A^, the time of core collapse is 
essentially independent of the initial binary fraction /, if / 10%, 
and smaller than the collapse time of a model without primordial bi- 
naries by a factor of roughly 2. This compares well with the expec- 
tation based on a two components model with siiigle stars only (and 
mass ratio 1:2) studied bv llnaeaki & Wivantd 1 1984), that would 
predict a relative ratio of the order 1.6 ± 0.2. While this lack of 
dependence on / may seem surprising, it should be observed that, 
at late times in core collapse, the binary fraction in the core has in- 
creased considerably (Fig. l2U . to 50% by number (or more) even if 
/ = 10% initially. Therefore the core is practically saturated with 
binaries, approximately independent of the initial value of /. Un- 
der these circumstances it is not surprising that the behavior of the 
core varies little with /. Further evidence of this "saturation" effect 
is visibl e below in F igs [l_8land ll9l 

Fre geau et alj J2003l) give a plot showing that the "core- 
collapse time" exhibits a roughly linear increase with /. This is 
quite at odds with our result, but the issue is a semantic one. 
iFreeeau et alJ |20o3) do not call the initial contraction "core col- 



lapse", preferring to reserve this term for the deep collapses akin to 
gravothermal oscillations; these occur much later, when the bina- 
ries are almost fully depleted. Thus the time they discuss might cor- 
respond to the deeper collapses in Fig.|5|after about t ~ ISOtr/i (0), 
as discussed in Sec. 14. II whereas our figure refers to the end of the 
initial phase of contraction at about t — 13trh(0). 



5.1.2 Core radius in the phase of steady binary burning 

The information from our simulations on this topic is summarised 
in Fig. llSI We notice again a saturation effect for initial binary frac- 
tions / 0.1, and a significant variation with A^. 

To explore these dependences further, we now compare the ra- 
tio of the core to half mass radius in the post-collapse phase with 
the theoretical estimate given by Vesoerini & Chernoff ( 1994). By 
balancing the energy production in the core with the rate of expan- 
sion at the half mass radius, they get the following estimate for this 
ratio (their eq. 9): 



4>b{i — 4>b)fj.bs + 4>bfJ-bb 



loB 



(3) 



,io(0.4iV) {l + (PbY 

where the various quantities are defined as follows. The quantity a 
is a parameter depending on 7 and Vc/v^, where 7 is the ratio of 
the expansion timescale rh/rh to trh, and Vc and Vh are the veloc- 
ity dispersion in the core and at the half-mass radius, respectively; 
we have assumed typ ical values for these parameters: 7 ~ 10 an d 
Vc/vh ~ ^/2 I Heggie & AarsetlJlT993 : lOoodman & Hutlll989l) . 
The quantity (pb is the binary fraction in the core defined as 

Nb 



Ns+Nb 



Finally, Hbs and Hbt are coefficients for the efficiency of binary- 
single and binary-binary burning and depend on the distribu- 
tion of binding energy of the binaries. These coefficients have 
been computed b y a numerical implementation of Equation 9 in 
IVesperini & Cher noff ( 1994), which has been checked against the 
results given by those authors. For simplicity, to compute the inter- 
action cross-sections we have assumed a flat distribution in log(i5f,) 
between 12 and 300 kT, where the upper and lower limits have 
been estimated from the measured binding energy distribution (see 
Fig.|9}. 

For a fixed initial binary fraction / = 0.1 this model is com- 
pared with the results of our simulation in Fig. llSl ta is assumed 
to be constant). Clearly the A'^-dependence is different. This result 
appears indeed rather puzzling. In fact it is unlikely that differences 
in the current fraction of binaries in the core, (j)b, are responsible, 
because of the saturation effect. The Hbs and fibb coefficients are 
also TV-independent, since we have checked that the distribution 
of binary binding energy, that is the only relevant input, evolves re- 
markably similarly in simulations with A'^ from 512 to 16384. In the 
worst case, our simplified approach (i.e. assuming a uniform distri- 
bution in log(_E6) for computing the analytical model) thus intro- 
duces only a constant multiplicative bias in eg. [3] In additio n, this 
bias should be limited, since Figs. 5 and 6 in lVesperini & Chernofj 
( 1994) suggest that the core radius is not particularly sensitive to 
the upper and lower limits of this distribution, at least within the 
limits of variation we observe (Fig. |9}. Other factors which may 
influence the efficiency of binary heating are the total mass and 
central potential, but M varies with A'' by at most about 2% up 
to t/trh{0) = 100 and (j>o does not show any systematic trend 
with A'^ (with an average value, after the core collapse, of « 1.6 in 
NB0DY6 units for our runs with / = 10%). 
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An interesting issue, raised by the referee, is tlie time scale 
on which binaries release energy. The time scale for close inter- 
actions between a binary of semi-major axis a and single sta rs of 
space density n may be estimated from the cross sections of iHutI 
il984) as approximately ths — Vi„f /{2AnGma), where Vinf is 
the relative velocity of a binary and third body when far apart, and 
m is the stellar mass. Comparing with the local relaxation time tr 
JSpitzer'l'igSf)'). we find that ^ 0.9Gm In A/(au^), where 

V is the root mean square three-dimensional velocity of single stars 
and In A is the Coulomb logarithm. Since this depends on A'^ there 
is an A'^-dependent regime of binaries which are effectively inert 
on the time-scale of a given simulation, and this potentially intro- 
duces an A^-dependence in the efficiency of binary interactions. On 
the other hand the time scale for interactions involving our hardest 
primordial binaries is approximately 240tr . Since the central relax- 
ation time is much smaller than the half-mass relaxation time, we 
conclude that all primordial binaries have enough time to interact 
within the time scale we have considered (up to 100t,.h(0)). 

We have also measured the value for 7 in the simulations, 
co nfirming, immediately af ter core collapse, the values reported 
by Heaeie & Aarseth (1992), i.e. 7 ~ 12, and observing basically 
no A^-dependence. In this regard we note that Gao et al. measured 
7 = 7. From our set of simulations this value appears to be un- 
expectedly low, even when account is taken of the fact that their 
model is isotropic (cf. Fig. 2 in Takahashi 1996). Though we are 
unable to explain the discrepancy between the values of 7 it does 
offer an alternative explanation of the fact that our values of rc/r^ 
appear to disagree with those of Gao et al (cf. our Sec.4. 1). Accord- 
ing to the theory of Vesperini & Chernoff a change in the value of 
7 leads to a change in the value of r^jru consistent quantitatively 
with what is found. 

Finally, the differences in the value of rc/rn at core collapse 
could be due to variations of the central velocity dispersion Vc, 
which could have a significant impact, as a oc {vc/vhf'- From 
our simulations we do not observe a significant A^-dependence 
for this quantity, wh i ch a ppears to be close to the values quoted 
bv iGoodman & Hull (1989), {vc/vhf « 2. To be sure a lim- 
ited A^-dependence is present, but this goes in the opposite di- 
rection with respect to the one required to obtain a better match 
in Fig. 1181 so that if this trend is taken into account the com- 
parison with the analytical estimate for rc becomes worse: for 
10% primordial binaries we measure {va/vhf = 1.88 ± 0.03 
with N = 256, {vc/vhf = 1.95 ± 0.02 for A^ = 1024 and 
[vc/vhf = 2.07 ± 0.02 for A^ = 4096. In contrast Gao et al. re- 
port a lower value, (vc/vhY ~ 1.8, so that this low central velocity 
dispersion could account for the smaller core in their Monte Carlo 
run, although it is not clear what is the origin of such a low velocity 
dispersion. Further studies would be required to clarify this point. 

To summarize, from the d ata obtained in our si mulations each 
term that contributes to the 1 Vesperini & Chemofll |I994) model 
cannot explain why rc/r^ exhibits a steeper variation in A' than 
expected. Since the number of particles for the lower A^ tail of the 
curve in Fig. ^| is rather modest, one possible explanation could 
be due to the importance of granularity effects so that a mean field 
approximation is not fully justified (see also the concluding discus- 
sion in Sec.|6j. 

At fixed A^ a comparison with eq.js} is shown in Fig. ll9l Here 
Q is kept as a free parameter. In this case, there is satisfactory quan- 
titative agreement within the scatter of different simulations at each 
value of /, except possibly for the highest values. (Note that we 
here use the current value of the binary fraction, <^j,, averaged over 
the same time interval, and not the initial value. The two groups 
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Figure 17. Time for core collapse for the set of simulations with 0% to 
50% primordial binaries. The filled tiiangles con'espond to / = 0. 



on the right, i.e. at ~ 0.8 and 0.9 correspond to initial val- 
ues of / = 0.5 and 1, respectively.) For A' = 2048 (not shown) 
we obtain comparable agreement for all /, except for an indication 
that the theory predicts slightly too small a core radius for the very 
highest and the very lowest values of /. It must be borne in mind, 
however, that the core radius in a multi-component system is not a 
well defined concept. Our data for A' = 2048 do not exhibit the 
slightly puzzling clumping at <f)b — 0.4 in Fig. 1191 or the system- 
atic offset at ~ 0.6. Our dataset for TV = 2048 is considerably 
larger, and we conclude that these features in Fig. ll9l are results of 
the more modest sample size. 



5.2 The binary and single populations 

5.2.1 Total mass in singles and binaries 

When the initial binary fraction is 10%, the binaries are mostly de- 
stroyed in the first lOOtrh(O). By the end of the simulation they 
account for only about 5% of the mass (Fig.0. The situation is 
quite different for a simulation which contains only binaries ini- 
tially (Fig. l20> . By the same point in this simulation the binaries 
account for about 80% of the remaining mass. 

5.2.2 Spatial Distribution of Binaries 

The concentration of binaries in the core, after the end of core 
collapse, is presented as a function of A' and / in Fig. 1211 Some 
comparison is possible with the results of Gao et al. at fixed A' 
(their Fig. 3b), who found that the ratio increased by 0.3 dex from 
/ = 0.1 to / = 0.2. An increase by this factor is consistent with 
what we find for all A' that we studied, but we repeat that the actual 
values we find are much smaller than those displayed by Gao et al. 
(Sec.E^J. 

While Fig. 12 II confirms that the concentration of binaries in 
the core increases with the original binary fraction, /, their overall 
spatial distribution differs less from that of the single stars, as we 
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Figure 18. Core to half mass radius averaged over a period of lOtrho ^fter 
core collapse and over different realizations (see TablelT] for the set of sim- 
ulations with 0% to 50% primordial binaries. The points associated with 
runs with primordial binaries are compared to the Vesperini & Chemoff 
model with 7 = 10, fc/fft = \/2 and (pf, = 0.4. The po ints associated 
with single-stars runs are compared to a N~^/^ power law iHeggie & Hull 
l2003h. Box 28.1 
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Figure 19. Core to half mass radius averaged for lOt^^o ^fter the core 
collapse for the set of simulations with A' = 4096 compared with the 
Vesperini & Chemoff model. For this model (Eq.|3| we have adopted 7 = 
9.4, tic /vfi = y/2, distiibution in the binary binding energy flat in log scale 
from 12 to SOOfcT. These values imply /.t^s = 2.841 and fi^b = 0.774. 



Figure 20. Dependence of number of single and binary stars (relative to 
the initial values) on time (units of the initial half-mass relaxation time). 
Simulation starting with 4096 primordial binaries and no singles. 
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Figure 21. Ratio of the number of binaries to singles in the core averaged 
for lOtr/iO ^fter the core collapse for the set of simulations with 10%, 20% 
& 50 % primordial binaries. 



consider models with larger /. This result is illustrated in Fig. 1221 
which should be compared with Fig. llOK se e caption to Fig.|22) . In 
this respect our result differs from that of Freseau et alJ as 
we have already remarked (Sec. 14.4k when / = 10% they find 
that the half-mass radius of binaries exceeds that of singles for 
t/trh{0) Si 120. This does not happen in their run with 20% bi- 
naries initially (their Fig. 5), where r^^t < fs^t at all times up to 
the end of their run, at t/trhiQ) — 600. 
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Figure 22. Dependence of half-mass radius (N-body units) for singles (up- 
per curve), whole system (second curve) and binaries (central curve) on 
time (units of the initial half-mass relaxation time). The lowest curve is the 
core radius of the system. The simulation has been performed with 4096 
particles and 50% of binaries. Unlike Fig. 1101 here, due to the higher ratio 
of primordial binaries, the half mass radius for binaries expands steadily 
after a short initial contraction. 

6 CONCLUSION 

This paper has had several objectives. Our first has been to pro- 
vide evidence from A'^-body simulations of the evolution of highly 
idealised stellar systems containing an initial population of bina- 
ries. The Fokker-Planck models of Gao et al. have become an im- 
portant touchstone for subsequent modelling, but until now it has 
not been clear whether any disagreement with their results is due 
to the simplifying assumptions they made by adopting an approx- 
imate method. What we have tried to do with our simulations is 
to establish what actually happens, so that subequent research can 
rest on fewer assumptions and approximations. In this respect our 
work can be considered as a stepping stone to improve the phys- 
ical understanding of the results from more complex and realis- 
tic numerical simulations, such as those recently performed by 
IPortegies Zwart & McMillanI j2004) . 

In line with this objective we have presented detailed results 
following closely the exposition of Gao et al., which covers the 
evolution of the structure of the system (core and half-mass radius, 
and structural parameters such as the scaled central potential), and 
the evolution of the binaries: their total numbers, their abundance 
in the core, their energies, and so on. 

Our second objective has been to begin a discussion of the 
differences between the various models which have appeared in 
the literature so far, with pa rtic ular reference t o those of G ao et 
aL. lGiersz & SpurzenJ <2000l) and lFregeau eTall J2003ll2005h . Here 
our discussion has covered several of the above issues, and also 
the gravothermal oscillations which have been observed in almost 
all simulations, but with somewhat contradictory results for some 
characteristics, such as the time of onset of the oscillations. We 
have also compared our results with some of the small amount of 
quantitative theoretical predictions in the literature, mainly those 



of IVesperini & Chemofj il994b . With respect to this analytical es- 
timate for the core to half mass radius in the post-collapse phase 
we find good quantitative agreement (at any fixed value of N) for 
the dependence on the initial binary ratio. On the other hand the 
observed A^-dependence of this quantity is much steeper than one 
would expect on the basis of the theoretical arguments given by 
f/esDerini & Chernoff 1 1994j). This could be simply a low-TV effect, 
but possibly some of the assumptions underlying the model could 
be inaccurate. 

This issue is relevant to the comparison of the core radius be- 
tween our models and Fokker-Planck or Monte Carlo models. With 
the present data we cannot give a clear answer to this point, since 
either extrapolation or direct simulations with TV > 10'' would be 
required. A different way of gaining some insight into this problem 
could be to run some comparison simulations with Monte Carlo 
methods with TV in the range considered in this paper. If our re- 
sults would be reproduced, they could be used to indirectly validate 
Monte Carlo simulations with TV high enough to answer this open 
question. 

Our third objective has been to attempt to add some new in- 
sight and information on the nature of the coUisional evolution of 
stellar systems with an initial population of binaries. In particu- 
lar we have considered the time of core collapse (the definition of 
which we have attempted to clarify). For the period of steady binary 
burning which follows this initial collapse, we have determined the 
manner in which the core properties (core radius and binary frac- 
tion) depend on the initial binary fraction and particle number. 

It is not clear which of these conclusions are applicable, even 
qualitatively, to real star clusters, with an initial mass function, a 
largely unknown initial distribution of binary para meters, a tidal 
field, and stellar and binary evolution. For example. Ilvanova et alJ 
( 2005) have argued (on the basis of a simplified treatment of the 
core) that the binary fraction in the core of a dense cluster may be 
as little as 5%, even if / = 1 initially. Now lFregeau et all j2003) 
find (with a Monte Carlo code) that a system in a tidal field with 
/ ^ 0.1 initially disrupts after at most 45frfe(0). If we were to 
guess from our Fig.|20| we would conclude that the binary fraction 
at this time would actually exceed 70%. Of course this argument 
ignores the essential physical issues of stellar evolution, collisions, 
etc., which underly the conclusion of Ivanova et al. (2005). On the 
other hand, all these models neglect factors which may be essential. 
Our next paper on this topic will provide detailed information on 
one of these which we have ignored so far - the influence of the 
tide. 
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